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THREE-DIMENSIONAL ELASTIC STRESS AND DISPLACEMENT ANALYSIS OF 
TENSILE FRACTURE SPECIMENS CONTAINING CRACKS 
by John P. Gyekenyesi, Alexander Mendelson, and Jon Krlng 

Lewis Research Center 

SUMMARY 

The line method of analysis is applied to the Navier- Cauchy equations of elastic 
equilibrium to calculate the displacement distributions in frequently used tensile frac- 
ture specimens. The application of this method to these equations leads to coupled sets 
of simultaneous ordinary differential equations whose solutions are obtained along sets 
of lines in a discretized region. When decoupling the equations and their boundary con- 
ditions is not possible, the use of a successive approximation procedure 'permits the 
analytical solution of the resulting ordinary differential equations. The results obtained 
show a considerable potential for using this method in the three-dimensional analysis of 
finite geometry solids and suggest a possible extension of this technique to nonlinear 
material behavior. 


INTRODUCTION 

Considerable progress has been made in recent years in the stress analysis of 
bodies containing flaws or cracks. However, most of the work on this subject has been 
based on the plane theory of elasticity (ref. 1). For a fracture specimen of finite geom- 
etry the stress and displacement fields are highly three-dimensional. Hence, solutions 
of the general elasticity field equations must be obtained when more reliable results are 
needed. 

Because of the geometric singularity associated with any crack type problem, there 
is almost no possibility of a simple closed form type of solution. For this reason, 
three-dimensional elastic solutions have been obtained only for a restricted class of 
problems (refs. 2 to 4). Recently, with the availability of large digital computers, a 
variety of numerical methods appeared in the literature; however, most of these meth- 
ods have yielded only partial results. Among these approximate methods are the finite 



difference (ref. 5), the direct potential (ref. 6), the eigenfunction expansion (ref. 7), 
and the line method of analysis (ref. 8). Of all these solution techniques, the line meth- 
od of analysis appears to yield the most complete and accurate results for three- 
dimensional elasticity problems. 

Although the concept of the line method for solving partial differential equations is 
not new (ref. 9), its application in the past has been limited to simple examples (ref. 10). 
The line method lies midway between completely analytical and discrete methods. The 
basis of this technique is the substitution of finite differences for the derivatives with 
respect to all the independent variables except one for which the derivatives are retain- 
ed. This approach replaces a given partial differential equation with a system of simul- 
taneous ordinary differential equations whose solutions can then be obtained in closed 
form. These equations describe the dependent variable along lines which are parallel 
to the coordinate in whose direction the derivatives were retained. Application of the 
line method is most useful when the resulting ordinary differential equations are linear 
and have constant coefficients. 

An inherent advantage of the line method over other numerical methods is that good 
results are obtained from the use of relatively coarse grids. This use of a coarse grid 
is permissible because parts of the solutions are obtained in terms of continuous func- 
tions. Additional accuracy in normal stress distributions is derived from the fact that 
they are expressed as first- order derivatives of the displacements and these derivatives 
can be analytically evaluated. Inherently inaccurate numerical differentiation is re- 
quired only for evaluating the shear stresses, but this presents no important loss of 
accuracy since they are an order of magnitude smaller than the normal stresses. For 
problems with geometric singularities, additional accuracy is derived from using a 
displacement formulation since the resulting deformations are not singular. 

It is the purpose of this report to present a simple and systematic approach to the 
elastic analysis of three-dimensional, finite geometry solids containing traction- free 
central or edge cracks. The need for these specific solutions has existed for a number 
of years in fracture toughness testing, and several attempts in the past have failed to 
arrive at meaningful results (refs. 11 and unpublished data by G. C. Sih and R. J. 
Hartranft of Lehigh University (NASA Grant NGR-39-007-025)). Problems that are most 
conveniently described in rectangular Cartesian coordinates are discussed herein while 
circular geometry solids will be treated in a subsequent report. 


SYMBOLS 

Ax 

A^ i = j =2, partitioned submatrices of e 

A(x) coefficient matrix of first-order, x-directional differential equations 
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2a 

B i 

2b 

E 

e 

G 

\> hyj \ 

I 

K I 

2L 

l 

m 

NX, NY, NZ 

n 

R 

r(x) 

s(y) 

2t 

t(z) 

u 

V 

w 

X, y,z 

x 


f] 


V V 1 

CT z’ °xy J V 
a zx’ CT yzJ 


crack width 

partitioned particular integrals, i = 1,2, 3, 4 
bar width 

Young's modulus of elasticity 
dilatation or natural logarithm base 
shear modulus of elasticity 
increments along Cartesian coordinate axes 
identity matrix 

stress intensity factor for opening mode 

coefficient matrices of second- order differential equations 

bar length 

number of lines in x- direction 
number of lines in y- direction 
number of lines in a given plane 
number of lines in z- direction 
distance from crack edge 

coupling vector for x- directional second- order differential equations 
coupling vector for y-directional second- order differential equations 
bar thickness 

coupling vector for z-directional second- order differential equations 

x-directional displacement 

y-directional displacement 

z-directional displacement 

rectangular Cartesian coordinates 

Lame's constant 

variable of integration 

Poisson's ratio 

components of the stress tensor in Cartesian coordinates 
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2 

V Laplacian in Cartesian coordinates 

Symbols used in the appendixes are defined there when they are introduced. 


REDUCTION OF THE NAVIER-CAUCHY EQUATIONS TO SYSTEMS 
OF ORDINARY DIFFERENTIAL EQUATIONS 

Within the framework of linearized elasticity theory, the equations of elastic equi- 
librium in terms of displacements are 


(x + G) — + GV 2 u = 0 
3x 


( 1 ) 


(X + G) — + GV 2 v = 0 

3y 


( 2 ) 


(x + G) — + GV 2 W = 0 
3z 


(3) 


where the body forces are assumed to be zero and the dilatation is 

e = + — 

3x 3y dz 

The stress- displacement relations can be written as 

E 


a x = 


CT„ = 


(1 + y)(l - 2v) 
E 


y (1 + !/)(l - 2v) 

— ! [( 

+ iv)(l - 2v) L 


[«■ 

[- 


a z = 


(1 


xy 


E /dv + 3u\ 
2(1 + u) \3x 3y ) 


(4) 


^)iH + „/iv + iwV] 

(5) 

3x \3y 3z/J 


V )& + v(*» + inV 

(6) 

3y \3z 3x/. 


„)iw + ,,/is + iiV 

(7) 

3z \3x 3y/ 



( 8 ) 
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( 9 ) 


a zx=-^— 
zx 2(1 + v) \3z 3x/ 

a = _ E _ /jw + 3v\ 
yZ 2(1 + v ) \3y 3z/ 


( 10 ) 


For a finite geometry solid with rectangular boundaries, we construct three sets of 
parallel lines (fig. 1(a)). Each set of lines is parallel to one of the coordinate axes and 
thus perpendicular to the corresponding coordinate plane. An approximate solution of 
equation (1) can then be obtained by developing solutions of ordinary differential equations 
along the x-directional lines. As seen in the figure, there are a total of l = NY x NZ 
such lines where NY is the number of lines along the y-direction and NZ is the number 
of lines along the x-direction in a given plane, respectively. We define the displacements 
along these lines as u^, Ug, . . . , u^. The derivatives of the y-directional displacements 
on these lines with respect to y are defined as v|j, v | • • •> v | ^, and the derivatives 

of the z-directional displacements with respect to z are defined as w L , w | g, . . . , w | 
These displacements and derivatives can then be regarded as functions of x only since 
they are variables on x-directional lines. When these definitions are used, the ordinary 
differential equation along a generic line ij (a double subscript is used here for sim- 
plicity of writing) in figure 1(b) may be written as 


A] , (1 - 2u) 


dx' 


2 2(1 - v ) 


(i + "ips + A 2 <<Ii+i -j + + h <Ui >i +i + "i.i-i’ 

VV h zJ h y h z 


f ij (x) 
2(1 - V ) 


= 0 


(H) 


where 


V x )=- 

1J dx 


ij 


dw 

dx 


ij 


v = 


dv 

dy 


(12) 


and 


w 


dw 

dz 
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Similar differential equations are obtained along the other x-directional lines. Since 
each equation has the terms of the displacements on the surrounding lines, these equa- 
tions constitute a system of ordinary differential equations for the displacements u^, 
u 2 , . . u r 

It is convenient to nondimensionalize equations (11) and (12) with respect to some 
characteristic dimension of the problem. Hence, we introduce the following variables: 





(13) 


where 2(a) is the crack width. In matrix notation, the differential equations along the 
x-directional lines of figure 1(a) can be written as 

d 2 u/dx 2 = u + r(x) (14) 

ixl ixi lx 1 lx 1 

In a similar manner, to solve equations (2) and (3) ordinary differential equations 
are constructed along the y- and z-directional lines, respectively. These equations are 
also expressed in an analogous form to equation (14); they are 

d 2 v/dy 2 = K y v + s(y) (15) 

mxl mxm mxl mxl 

d 2 w/dz 2 = K w + t(z) (16) 

’ Li 

nxl nXn nxl nxl 

An inspection of these equations shows that they are linear, second- order, ordinary 
differential equations with constant coefficients. The coupling terms in each set of dif- 
ferential equations are grouped into the vectors r(x), s(y), and t(z), which are the 
nonhomogeneous terms in equations (14), (15), and (16). The elements of the coefficient 
matrices are functions of the coordinate increments and Poisson's ratio only. Since the 
sum of the elements in any given row of the three coefficient matrices is zero, they are 
all singular. 
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Noting that a second- order differential equation can satisfy only a total of two 
boundary conditions and since three-dimensional elasticity problems have three bound- 
ary conditions at every point of the bounding surface, some of the boundary data must 
be incorporated into the surface line differential equations. Hence, conditions of normal 
stress and displacement are enforced through the constants of the homogeneous solutions 
while shear stress boundary data must be incorporated into the differential equations of 
the surface lines. The application of the specified shear conditions permits the use of 
central difference approximations when surface line differential equations are construct- 
ed. The details of constructing these equations are found in reference 8. 


SOLUTION OF THE SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 


Appendix A shows the details of a solution technique for the type of simultaneous 
ordinary differential equations given in equation (14). This solution technique yields the 
following solution vector for the x-directional displacements: 


u(x) = A n (x) 

u(0) + Bj(x) 

+ a 12 (x) 

u(0) + B 2 (x) 

X 

'"O 

H 

X 

1 

i-H 

X 

r>0 

X 

ixl 

jxl Zxl_ 


where A 1]L (x) and 


Aj^Cx) are the matrix functions 


I 

/ 

'i 


A 11 (x) = cosh K x / x 


(17) 


(18) 


a 1 2 (x) = K^ 1/2 sinh K^ /2 x (19) 

1/2 ~ 

The matrix is defined as a matrix whose square is equal to 1^. Vectors u(0) 

and u(0) are initial value vectors whose evaluation from given boundary conditions is 
discussed in appendix B. Vectors B 1 (x) and B 2 (x)> which represent particular solu- 
tions of equations (14), are given by 


~ r x 

Bj(x) = - / A l2 (T?)r(7]) d t] 

lx 1 JO ixl lx 1 


B 2 (x) = / A 1 j(rj)r(77) dy 

lx 1 Jo ixl lx 1 


( 20 ) 


( 21 ) 
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Differentiating equation (17) with respect to x yields 


/VI ,/>-> 

L 

J 



J 

u(x) = A 12 (x) 

u(0) + Bj(x) 

+ A n (x) 

u(0) + B 2 (x) 

lx 1 txi ixl 

[lx 1 

Zxl J 

ixl 

jxl 

JX1J 


(22) 


Similar solution vectors can be constructed for equations (15) and (16) involving the y- 
and z-directional displacements, respectively. Since rfa) in equations (20) and (21) is 
unknown, the solution of the problem is begun by assuming values for dv/dx, dv/dx, 
dw/dx, and dw/dx along the x-directional lines. Using zero values for these required 
quantities is a good place to start. Then equations (17) and (22) give the first estimate 


for the vectors u(x)^ and It is assumed, of course, that the boundary con- 

rv rsJ / \ /-**>/ 1 j 

ditions u(0) and u(0) are specified or known. Using these calculated values of u v ' 

~( 1 ) ; 
and u v ' we can evaluate the vector s(y) v ' where we must use similarly assumed 

values of dw/dy and dw/dy. Analogous equations to equations (17) and (22) give the 

~ /V. /I \ Aj />»» / 1 \ 

first value of v(y) v ’ and v(y) v First values of w(z) v ; and w(z) v ' can then be 
calculated by using the first estimates of the x- and y-directional displacements and 


id) 


their derivatives in the coupling vector t(z)d). At this point we return to equations (17) 


( 2 ) 


and u(x)d) based on 'the first values 


i(i), 


and (22) and calculate the second values of u(x)' 
of the y- and z-directional solutions. 

If the values of u(x)^, u(x)^), v(y)^) ? y(y)^\ w( z)^), and w(z)''*' converge 
with the repetition of this procedure, an approximate solution of a given problem is 
determined. In general, the convergence rate for these successive calculations is de- 
pendent on the accuracy to which the matrix functions can be evaluated. Equa- 
tion (A13) may be used to establish this required accuracy. Sufficiently large errors in 
the necessary matrix functions lead to divergence in the successive approximation cal- 
culations. 

Since the vector r (17) in equations (20) and (21) involves displacements and their 
derivatives that are defined only at the nodes, finite difference calculus must be used in 
evaluating its elements. Hence, integrals B^(x) and Bg(x) are calculated by a suitable 
numerical integration technique . 

Once the displacement field in the bar has been calculated and the successive ap- 
proximation procedure has converged, the normal stress distributions along the three 
sets of parallel lines can be obtained from the following equations: 


„ + — — 

<l + ,)U-2„) ^ es (l + ,)(l-2,) 

£Xl £Xl 


(v + w), 


along 

x-lines 


Jxl £xl 


(23) 
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/V /N/ 


a = — ~ — — (v) , + — (u + w) , 

y TSL ISfes 


(24) 


mxi 


mxi 


mxl mxi 


a - - E(Lt.. 1_ , + gg 

2 (1 + rt(l - 2 rt |i°gf es a + ■')<! - 2 y) 


nxl 


nXl 


(u + v) alon g 

x-lines 
nXl nxl 


(25) 


Note that the normal stress equations involve only derivatives that can be evaluated 
analytically. The shear stresses at each node can be calculated from equations (8) 
to (10) but finite difference approximations must be used for the required displacement 
gradients. 


APPLICATION TO TENSILE FRACTURE SPECIMENS CONTAINING CRACKS 

A great amount of experimental work has been done in fracture mechanics (ref. 12) 
through the use of crack- notched specimjens. In the past, many different types of speci- 
mens have been used to [determine a material's fracture toughness. The most common 
early specimens employed in these tests' were the center- cracked and double-edge- 
notcljed bar specimens. Figures 2(a) arid 7(a) show the finite rectangular bars with 
through- thickness, traction-free central and double-edge cracks, respectively. Be- 
cause of the symmetric geometry and loading, only one-eighth of the bars has to be dis- 
cretized as shown in figures 2(b) and 7(b). Inspection of the derived ordinary differen- 
tial equations also shows that for any numerical computation a value of Poisson's ratio 
must be selected. For our work, a Poisson's ratio of 1/3 is used as shown in figures 2 
to 14. 


NUMERICAL RESULTS 

Center-Cracked Tensile Fracture Specimen 

The solution of this problem was obtained by using two different sets of lines along 
the coordinate axes so that the convergence of the finite difference approximations could 
be checked. In a given direction, uniform line spacing was used in all computations with 
no other restriction being placed on the selection of the grid size. The crack edge loca- 
tion with respect to the imposed grid was assumed to be halfway between nodes specify- 
ing normal stress and displacement boundary conditions, respectively. The successive 
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approximation procedure required for decoupling the three sets of ordinary differential 
equations was terminated when the difference between successively calculated displace- 

_ c 

ments at every point was less than a preset value (10” ). 

Selected results of the dimensionless displacements are listed in tables I to ni. 

The dimensionless crack opening displacement is shown in figure 3. Inspection of fig- 
ures 3(a) and (b) shows that the grids of lines have been sufficiently refined when the 
results of the 48 by 96 by 128 grid were calculated. Figure 3(a) also contains the re- 
sults of the plane elasticity solutions obtained by Mendelson, Gross, and Srawley 
(ref. 13). It is noteworthy that the results correspond to elliptical crack profiles in all 
cases. The plane stress solution gives the highest crack opening displacement while the 
plane strain solution is very close to the results obtained at the center of the bar. 

The dimensionless normal stress distributions in the crack plane are shown in fig- 
ures 4 to 6 as a function of both the x- and z- coordinates. The results in these figures 
clearly indicate the singular nature of the normal stresses near the crack edge. As ex- 
pected, the stress normal to the crack plane increases most rapidly near the crack 
front. A plot of these stresses in the z-direction also indicates a central region of uni- 
form stress and a boundary layer through which the stresses decrease to the surface 
values. These same results show that as x increases the stress field approaches a 
uniaxial state of stress which indicates that the cause of a triaxial stres's field in the bar 
is the through-thickness central crack. 


Double-Edge Crack Tensile Fracture Specimen 

Selected results for the problem of figure 7(a) are presented in tables IV to VI and 
figures 8 to 11. These results were also obtained with two different sets of lines along 
the coordinate axes which were identical to the two sets used for the central crack prob- 
lem. The dimensionless crack opening displacements are plotted in figure 8. Com- 
paring figures 8(a) and (b) shows that the finite difference approximations have suffi- 
ciently converged when the results of the finer grid were calculated. Contrary to the 
central crack problem, the crack opening displacements in figure 8 are independent of 
the z- coordinate. Similar results were obtained by Cruse and Van Buren (ref. 6) for 
the single- edge- crack bar specimen. 

The dimensionless normal stress distributions in the crack plane, are shown in fig- 
ures 9 to 11 as a function of both the x- and z- coordinates. Their singular nature near 
the crack edge is evident from the distributions shown in these figures. A plot of these 
stresses in the thickness direction also shows a central region of uniform stress and a 
boundary layer through which the stresses decrease to the surface values. Tables IV 
to VI contain selected results of the dimensionless displacement distributions. 
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STRESS INTENSITY FACTOR 


It is customary in fracture mechanics to describe the plane elasticity crack opening 
displacement as a superposition of three basic deformation modes (ref. 1). Since the 
problems shown in figures 2(a) and 7(a) have geometric symmetry and are symmetri- 
cally loaded, only the opening mode of crack displacement is obtained. In terms of the 
stress intensity factor for the opening mode Kj, the plane elasticity crack displace- 
ments near the crack tip are given by (ref. 1) 

vL_ n = ^ K t \ — plane strain (26) 

y ~ u G 1 1 2n 


v 


2 

y=° ' a + v)q 


K t 



plane stress 


(27) 


Since three-dimensional problems are neither in a state of plane strain nor in a state of 
plane stress, the definition of a stress intensity factor for these problems must be first 
established. Note that by definition the plane stress and plane strain stress intensity 
factors are equal while the displacements are approximately 12. 5 percent different for 
v = 1/3. Since the results indicate that most of the bar in the thickness direction is 
approximately in a state of plane strain, equation (26) is selected to calculate the stress 
intensity factor. Rearranging this equation so that the dimensionless crack opening dis- 
placements can be used leads to 



(28) 


where 


Cj = - 4(1 (29) 

E^2na. 


A plot of equation (28) as VR/a — 0 can then be used to calculate Kj. Since the crack 
opening displacement is a function of the thickness variable, the previous stress inten- 
sity factor varies in the z-direction. However, if we were to account for the nonplane 
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strain conditions near the surface by using equation (27) or a corrected equation (26) for 
the definition of Kj, the stress intensity factor would become a constant across the 
thickness of the bar by definition. This approach would then result in a continuously 
varying definition of Kj across the thickness . It should be noted that the previous 
description of Kj is completely arbitrary and that it is questionable if it has any real 
significance in three-dimensional elasticity problems. However, values of Kj are still 
presented so that a comparison is possible between the calculated results and the pub- 
lished plane strain solutions (ref. 1). 

Figures 12 and 14 show the calculation of the opening mode stress intensity factors 
from the plane strain crack opening displacements. Their variation across the thick- 
ness of the bar is shown in figures 13 and 15 for the two problems discussed in this re- 
port. Note that for the center- cracked bar Kj is maximum at the surface and minimum 
near the center. The value of Kj for the double- edge- crack specimen is independent of 
the z- coordinate since the crack opening displacement is constant across the thickness. 

Although the stress intensity factors for these problems could be determined with 
reasonable accuracy, the associated type of singularities are difficult to evaluate be- 
cause values of the normal stresses are needed within a distance of 0. 05a or less from 
the crack edge. With the equal spacing of lines used in these examples, the minimum 
node location for these problems is about 0. 06a. For this range of crack edge distance 
R, the singularity of the stresses is not defined. 

Lewis Research Center, 

National Aeronautics and Space Administration, 
and 

U.S. Army Air Mobility R&D Laboratory, 

Cleveland, Ohio, January 4, 1973, 

501-21. 
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APPENDIX A 


SOLUTION OF THE SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 
WITH CONSTANT COEFFICIENTS 

For the solution of equations (14) the following new variables are introduced: 



In terms of these new variables, equations (14) can be written as a set of twice as many 
first-order differential equations. Hence, we have 

4 U = A U + r(x) (A2) 

dx 2ZX1 21X21 2lx\ 2Vx\ 

where 

”0 i" 
ixi ixi 

A = (A3) 

21X21 K x 0 
1X1 ixi 



The solution of equation (A2) is well known and can be written as (ref. 14) 


13 



/V 


U(x) = e A * uC6) + e A * 
2Zxl 2ZX2Z 2ixl 2Zx2£ 


f 


e~^ r (?]) d?] 
2ZX2Z 2ZX1 


(A5) 


where U(0) is an initial value vector whose evaluation from given boundary conditions is 

Z-'-' 

A x 

discussed in appendix B and e is a matrix series given by the following equation: 


a~ a~ » 2~2 a 3~3 

e Ax = i + Ax + Aj^ + Aj L + 


2 ! 


3! 


(A6) 


In terms of the coefficient matrix equation (A6) yields 



OO 


oo 



x 2111 * 1 ^->1 V x 2m r m 

(2m + 1)! ^ /_ u (2m)! ^ 


r i' 

ii 

o 


o 

II 

a 

e Ax = 

A n (x) 

ixl 

a 1 2 (x) 

ixi 

21x21 

A 21 (x) 

A 22 (x) 


I'Xl 

I'Xl 


From equations (A7) and (A 8) we note that 




Ajj(x) = A 22 (x) = cosh 


1/2 


x 


(A7) 


(A8) 


(A9) 


A 12 (x) = K^ 1//2 sinh K*/ 2 x (A10) 

a 2 i(x) = k x a i 2 (x) (All) 

where K^ 2 is a matrix whose square is equal to K x< Noting that A^ and A 22 are 
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even functions of x while and A 2 j are odd functions of x, the inverse of e Ax 
becomes 


A n (x) -A 12 (x) 

-A 2 i(x) A 22 (x)^ 

Substituting equations (A8) and (A12) into the identity of e Ax ■ e” x = I yields 

A ll ( x)- I^A^Cx) = I 
1*1 1*1 1*1 1*1 



(A12) 


(A13) 


Matrix functions (A9) and (A10) may be evaluated for each value of the independent 
variable from the series definitions given in equation (A7). However, for increasing 
values of x serious convergence difficulties may arise and unpractically large number 
of terms must be calculated. To avoid these computational difficulties, additive formu- 
las for these matrix functions may be obtained from using the identity 

Ax^ Ax 2 A^Xj+Xg) 
e • e = e 

In terms of the submatrices A^. this identity yields the following two equations: 

A ll^ x l + x 2^ = A 11^ x 1^ A 11^ x 2^ + ^x A 12^ x 1^ A 12^ x 2^ (A14) 

and 

A 12^ x l + x 2^ = A 12 ^ x 1^ A 11^ x 2^ + A 11^ X 1^ A 12^ X 2^ (A15) 

Using equations (A7) in conjunction with equations (A14) and (A15) does not require 
the solution of an eigenvalue problem. However, greater accuracy is obtained if we 
consider the similarity transformation of 


A =P _1 K x P (A16) 

1*1 1*1 1*1 1*1 

where A is a diagonal matrix whose elements are the eigenvalues of 1^. Appendix C 
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gives the details of the development showing how the eigenvalues and eigenvectors of a 
matrix having the form of can be analytically evaluated. The eigenvalues of 
are given by 


x ij " 2k 3 


1 - cos 


(— )"1 

\NZ - 1/ J 


+ 2kr 


1 - cos 


j - 1 
NY - 1 


(A17) 


where i = 1,2, . . NZ, j =1,2, . . ., NY, and kg and kg are constants given in 
appendix C. Note that these eigenvalues are ordered by fixing i first and then varying 
j- 

The modal matrix or the matrix of eigenvectors corresponding to the previously 
ordered eigenvalues is given by 


P - P 2 ® Pj (A18) 

ixi NZXNZ NYXNY 

where ® denotes the Kroenecker product of two component modal matrices (ref. 16) 
whose elements are given by 


p ! 

= cos Dii.r 

— * s, j =1,2, . 

• , NY 

(A19) 


NY - 1 



P 2 

= cos < r - D# - 

— ^ r , i = 1,2, . 

. . , NZ 

(A20) 

^ri 

NZ - 1 


1 


The matrix functions (A9) and (A10) can now be evaluated in diagonalized form and re- 
transformed according to the transformation 


ii * PA n p 1 

(A21) 

12 “ PA 12 P 1 

(A22) 


Note that the inverse of the modal matrix P can also be evaluated in closed form and 
the details of the derivations are shown in appendix C. 

In the explicit evaluation of equation (A5), the value of the particular integral cannot 
be obtained until the vector r(Tj) is known along the x- directional lines. We define 
column vectors B^(x) and Bg(x) as follows: 
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(A23) 


|b 2 (x)J 
2Zxl 

In terms of the partitioned matrix functions, the integral in equation (A23) becomes 


~ f x 

Bj(x) J Q A 12 (T?)r()j) drj 

(A24) 

~ f* 

B 2 (x) = J A n (r])r(r]) dr] 

(A25) 


'0 


/v 

-X 


f e~^ r ( 77 ) drj 


21*21 21*1 


If it is assumed that r(r]) is known, the particular integrals (eqs. (A24) and (A25)) can 
be evaluated. Using the partitioned form of the matrices we can write equation (A5) as 


p(x) 



r 


A n (x) A 12 (x) 
A 2 i( x ) A 22^ x 1 


u(0M jBjtxjj 

u(0)[ ]b 2 (x)J 


(A26) 


from' which the solution vectors (17) and (22) can be directly constructed. 
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APPENDIX B 


DEVELOPMENT OF INITIAL VALUE VECTORS 

Since the problem in figure 2(a) is a two point boundary value problem, the initial 
values of both u and u are usually not available. A method for evaluating the initial 
value vectors for equations (14) and (15) from given boundary conditions is now pre- 
sented. 

From symmetry conditions, we can immediately conclude that 


Fj_ = u(0) = 0 (Bl) 

ixl zxl Zxl 


The zero normal stress boundary condition on the face x = b is used to evaluate the 
vector u(0). From equation (22) we have at x = b 



Using equation (5) with the given normal stress boundary condition gives 


(B2) 



- v 

1 - V 



(B3) 


Provided A 1;l (b) is not singular and using equations (Bl) and (B3) we find from equa- 
tion (B2) that 

F 2 = u(0) = -^- A‘}(b)(v + w) - AjJ(b)A 21 (b)B 1 (b) - B 2 (b) (B4) 

Zxl Zxl ixl lx 1 ixi ixl zxl zxl 


Similar equations are obtained for the initial value vectors of equations (16) along the 
z-directional lines. 

Along the y- directional lines the boundary conditions are somewhat more involved 
since in the crack plane mixed boundary conditions are specified. Denoting the number 
of y- directional lines falling over the crack surface as NIC and those falling outside as 
NOC, the zero normal stress condition over the crack face and the symmetry condition 
in the crack plane result in the following equations: 


/■«j 

v(0) 

NICXl 


- v 




(u + w)~ 0 


1 - V 

inside crack 


(B5) 
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v(0) = (0)„ (B 

NOCxl y_0 

outside crack 

If we assume that on the face y = L there is a uniform stress of ct q as shown in fig- 
ure 2(a), the normal stress boundary condition on this face gives 


v( L ) .fg a + ^arj-) . ( ~ + j 

E (1 - v) 1 - v y=L 

mxl mxl mxl mxl 


This vector can be suitably partitioned into vectors v a (L) and v^(L). For convenience 

NICxl NOCxl 

of matrix manipulations, the following vectors are constructed: 


3a 

NICxl 


v(0) 

mxl 

Moil J . 

1 J v(0) 
2mxl m xl 


F 

o/3 

NOCxl 


4a 

NICxl 


r 4/3 
NOCxl 


For the central- crack problem, values of F^ a and Fg^ are given by equations (B5) 
and (B6), respectively. An analogous solution to equation (A26) along the y-directional 
lines can be written as 


v(y)| [b n (y) D i 2 (y)| /p(°)| fB 3 (y) 


v(y) D ?1 (y) D 99 (y) \\{ v(0) B 4 (y) 


From equation (B9), we can express v(L) in a partitioned form consistent with equa- 
tion (B8) as 
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r- 

- 


r 


r 

\ 

v a 


D 21al 

D 21a2 

/ 

F 3a 


B 3* 

\ 

NICxl 


NICxNIC 

NICXNOC 


NICxl 


NICxl 


< 

> 



■« 


> + < 


> 

% 


D 2 1/31 

D 21/32 


F 30 


B 3/3 


NOCxl 


NOCXNIC 

NOCXNOC 


NOCxl 


NOCxl 



y=L 



y=L 

\S J 



y=L / 


+ 


°22al 

D 22a2 

/ 

( 

a 

Cn 


f 

B 4q? 

NICXNIC 

NICXNOC 


NICxl 


NICxl 




< 

► + < 

> 

D 22/31 

D 22/32 


F 4/3 


B 4/3 

NOCXNIC 

NOCXNOC 

y=L 

NOCxl 

V 1 - J 


NOCxl 

v_ J 


J 


(BIO) 


Equation (BIO) leads to two matrix equations involving the two unknown vectors 
and F^. The solution of these equations yields 


4/3 


>L‘ 


/3 


y=L 


D a D 21/31 


D 


y=L 


2l0!l 


y=L 


a 


y=L 


NOCxl NOCXNOC NOCxl NOCxNOC NOCXNIC NICxNIC NICxl 


- B 


4/3 


d; 


-l 


y=L 


D b B 3/3 


y=L 


NOCxl NOCxNOC NOCxNOC NOCxl 


,-l 


D_ 


' F 4q; + B 4a 


y=Ly 


NOCXNOC NOCxNIC NICxl NICxl 


(Bll) 


where 


D„ 


NOCXNOC 


( D 22/32 " D 21/31 D 21q- 1 D 22a2) 

y=^ 

NOCXNOC NOCXNIC NICxNIC NICXNOC 


(B12) 
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D b " ( D 21/32 " D 21/31 D 21ofl D 21a2) 

y=L 

NOCxNOC NOCXNOC NOCXNIC NICxNIC NICxNOC 


(B13) 


D c " ( D 22/31 ” D 21 j 31 D 2lQ!l D 22al), 

NOCXNIC NOCXNIC NOCXNIC NICxNIC NICxNIC 


y=L 


(B14) 


F 3a - D 21o!l 

/■v> 

„ ~ ' D 2lQ!l 

~ ^ D 2lQ!2 

~ ~ B 3/3 

B 3a 


y-L 

y=L 

y=L 

y=L 

y=L 


y=L 


NICXI NICXNIC NICX1 NICxNIC NICxNOC NOCxi NICxl 


D 21q-1 

„ ~ D 22al 

„ „ ( F 4a + B 4a 


y=L 

y=L \ 


y=Ly 


NICxNIC NICXNIC NICxl NICxl 
.-1 


'21 al 


y-L 


22ff2 


y-L 


F 4)3 + B 4j3 


y=L> 


NICXNIC NICXNOC NOCxi NOCxi 

(B15) 


with given by equation (Bll). Note that although the full matrix Dgj 


is 


r*j 

y=L 


singular, the partitioned submatrices are not. Equations (Bll) and (B15) together with 
the given boundary data completely specify the initial value vector needed for the solu- 
tion of equations (15). 
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APPENDIX C 


EVALUATION OF THE COEFFICIENT MATRIX EIGENVALUES AND EIGENVECTORS 

For the x-directional lines shown in figure 1(a), the nondimensionalized coefficient 
matrix K x can be arranged as 


NYXNY NYXNY 


K x 2 K xl K x2 
NY<NY NYXNY NYXNY 


^ = 


x2 *xl x2 

NYXNY NYXNY NYXNY 

0 2K x2 K xl 

NYXNY NYXNY 


where the submatrices K xl , K x g, and their elements kj, kg, and kg are 

k l = 2 ( k 2 + k 3 } 1 


k _ (1- 2u) 1 

2 2(1 - v) 1^2 

\ y/ 


2(1 - u)\r 2 
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A close investigation of equation (Cl) shows that the coefficient matrix K^. can be 
decomposed into component matrices having the following tridiagonal format: 



MxM 

The eigenvalues and eigenvectors of this type of matrix can then be obtained in closed 
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form. Noting that in equation (Cl) we have NZ rows of submatrices each of order NY, 
we express the coefficient matrix as 

= k 2 I 1 ® Kj + k 3 K 2 (x) I 2 (C6) 

ZXZ NZXNZ NYXNY NZXNZ NYXNY 

where (x) denotes the Kroenecker product of two matrices (ref. 15). Matrices Kj and 
K 2 have the desired form of matrix (C5) but are of different order. Associated with the 
matrices Kj and K 2 are the following two eigenvalue problems: 

XI = ju XI (C7) 

NYXNY NYX1 NYX1 

K 2 X2 = 6 X2 (C8) 

NZXNZ NZxl NZX1 

where jUj> j = 1,2, . . .,NY denote the eigenvalues of Kj and 6j, i = 1,2, . . .,NZ 
represent the eigenvalues of K 2 - The original eigenvalue problem associated with the 
coefficient matrix can be written as 

X3 = AX3 (C9) 

ZXZ lx 1 zxl 

After some matrix manipulations involving Kroenecker products (ref. 15), it can be 
shown that the eigenvalues A- and the corresponding matrix of eigenvectors X3^ can 
be expressed in terms of the component matrix eigenvalues and eigenvectors. The re- 
sults of these manipulations are 


X ij ' k 3 5 i + Vj 

(CIO) 

X3^ = X2 1 ® Xl^ 

(Cll) 

ZXZ NZXNZ NYXNY 



where i = 1, 2, . . .,NZ and j = 1, 2, . . .,NY. Equations (CIO) and (Cll) reduce the 
problem of equation (C9) to that of finding the eigenvalues and eigenvectors of ma- 
trix (C5). 

The eigenvalues of the tridiagonal matrix (C5) can be obtained by using difference 
equation theory. The details of this method for finding the eigenvalues of matrix (C5) 
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can be found in reference 8, but the results for the problem of equation (C7) are 



2 


1 - cos 



7 T 


j = 1,2, . . . , M 
M = NY 


(C12) 


Similarly, the eigenvalues of Kg become 


1 - cos 

VI 

i = 1,2, . 

. . , M 

(C13) 

\M - 1 

■/ J 

M = NZ 




Substituting equations (C12) and (C13) into equation (CIO) leads to the results of equa- 
tion (A17). 

The details of finding the eigenvectors corresponding to equations (C12) and (C13) 
are also given in reference 8. The results of these manipulations are given by equa- 
tions (A19) and (A20), respectively. Substituting these eigenvectors into equation (Cll) 
leads to equation (A18) which gives the modal matrix of the coefficient matrix K x . 

The inverse of the modal matrix P can be obtained from 

P _1 = Pg 1 0Pj 1 (C14) 

provided that closed form inverses of P-^ and Pg can be found. Since these two ma- 
trices are of the same type, it is sufficient to investigate the closed form inverse of Pj 
only. 

As the first step in developing the closed form inverse of Pp we construct a diag- 
onal matrix which transforms the nonsymmetric matrix Kj into a symmetric 
matrix K®. The form of this transformation is 

d 1 /2k 1 d ; 1/2 = K® (C15) 

where by definition the square root of a diagonal matrix is also a diagonal matrix whose 
elements are the square roots of the elements in the original matrix. For the tridiag- 
onal matrix Kj of order M, the diagonal matrix Dj is of the form 
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D 1 = 
MXM 


1/2 0 0 0 0 0 O' 


0 0 0 0 


0 0 1 0 0 0 0 


0 0 


0 0 0 0 0 


0 0 


0 0 0 0 0 0 1/2 


Associated with matrix is the similarity transformation of 


(C16) 


K 1 P 1 = P 1 A 1 (C17) 

which can be written as 

KjDJ^dJ^Pj = P 1 A 1 (C18) 

1/2 

Premultiplying equation (Cl 8) by Dj' gives 

K 1 ( D l /2p i) = ( D l /2p i) A 1 (C19) 

Since K® is a symmetric matrix, the transformation of equation (C19) is orthogonal and 
the columns of (dJ^Pj) must be orthogonal. For an orthogonal transformation it is 
known that 



(C20) 


where T denotes the transpose and D is a diagonal matrix. Since D x is a diagonal 
matrix, the first term in equation (C20) can be written as 



- P T D 1/2 

- ^ u x 


(C21) 
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Using equation (C21) in equation (C20) yields 


P 1 D 1 P 1 =D 


Premultiplying this equation by D~ * gives 



(C22) 


(C23) 


from which we find that 


Pj 1 =D _1 PjD 1 (C24) 

Equation (C24) requires the inverse of an unknown matrix D. However, from equa- 
tion (C22) it can be shown that the matrix D always has the form 


D = 
MXM 


M - 1 0 0 


0 M.._ 1 0 0 0 


0 0 M " 1 0 


0 0 0 0 M - 1 


Comparing matrices (C25) and (C16) we can conclude that. 


d = M ~ 1 d: 1 

O «*■ 


(C25) 


(C26) 


Substituting this equation into equation (C24) gives the final closed form inverse of Pj 
as 
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2 


(C27) 



MXM 


M - 1 



MXM MxM MXM 


A similar equation can be written for the inverse of the modal matrix Pg. With the 
closed form solution of the component modal matrices and their inverses, the diagonal- 
ized form of K x is obtained from equation (A16). 
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T z 70 by 98 grid), 

lb) Discretized region of rectangular bar with through -thickness central crack. Figure 3 . . Dimension|ess crack opening djsp|acement 

for rectangular bar under uniform tension containing 

Figure 2. - Rectangular bar with through-thickness central crack under uniform tension. through-thickness central crack. 
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Figure 4. - Dimensionless x-directional normal stress distribution Figure 5, - Dimensionless y-directional normal stress distribution 

in crack plane for rectangular bar under uniform tension con- in crack plane for rectangular bar under uniform tension con- 
taining through-thickness central crack. taining through-thickness central crack. 



.'Crack edge location 
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a =1.0 v = 1/3 
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ro o 
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(a) Dimensionless crack opening displacement (48 by 
96 by 128 grid). 



(b) Dimensionless crack opening displacement (35 by 70 
by 98 grid). 

Figure 8. - Dimensionless crack opening displacement 
for rectangular bar under uniform tension containing 
through-thickness double-edge cracks. 



(a) Dimensionless x-directional normal stress as 
function of x. 



(b) Dimensionless x-directional normal stress 
as function of z. 

Figure 9. - Dimensionless x-directional normal stress 
distribution in crack plane for rectangular bar under 
uniform tension containing through-thickness 
double-edge cracks. 
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a - 2.0 
b = 2.0 
L = 1. 75 


t- 1.5 
v- 1/3 



Figure 12. - Calculation of stress intensity factors Kj 
for rectangular bar under uniform tension containing 
through-thickness central crack. 
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Figure 13. - Variation of stress intensity 
factor Kj across thickness for rec- 
tangular bar under uniform tension 
containing through -thickness central 
crack. 



R/a 

Figure 14. - Calculation of stress intensity factors Kj 
for rectangular bar under uniform tension containing 
through-thickness double-edge cracks. 
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Figure 15. - Variation of stress intensity 
factor Kt across thickness for rec- 
tangularlwr under uniform tension 
containing through-thickness double- 
edge cracks. 
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